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SYSTEM, METHOD AND APPARATUS FOR SMALL PULMONARY 
NODULE COMPUTER ATOED DIAGNOSIS FROM 
COMPUTED TOMOGRAPHY SCANS 

This application claims the benefit of U.S. Provisional Application 
No.60/41 9,597, filed October 18, 2002, which is incorporated herein by reference. 



BACKGROUND OF THE INVENTION 

5 The present invention relates to the art of diagnostic imaging of small 

pulmonary nodules. In particular, the present invention is related to analyzing and 
manipulating computed tomography scans to: detect lung nodules and correlate a pair 
of segmented images of a lung nodule obtained at different times. 

10 Lung cancer is the leading cause of cancer death in the United States. 

According to the American Cancer Society, there were approximately 169,400 new 
cases of lung cancer (90,200 among men and 79,200 among women) in the United 
States in the year 2002. About 1 54,900 lung cancer deaths were predicted for the 
same year, which accounts for 28% of all cancer deaths. Although survival rate of 

1 5 lung cancer is only 14%, results from the ELCAP project show that detection and 
treatment of lung cancer at early stages may improve this rate to 70%. C. I. 
Henschke, D. I: McCauley, D. F. Yankelevitz, D. P. Naidich, G. McGuinness, O. S. 
Miettinen, D. M. Libby, M. W. Pasmantier, J. Koizumi, N. K. Altorki, and J. P. 
Smith, "Early Lung Cancer Action Project: overall design and findings from 

20 baseline screening." Lancet July 10, 1999; 354(9173):99-105. Hence, lung cancer 
screening has recently received considerable attention. 



25 



In the screening process, radiologist analyze images of asymptomatic 
patients searching for a specific abnormality. Henschke et al reported that using 
low-dose CT as compared to chest radiography can improve the detection of small, 



non-calcified nodules at potentially more curable stage. Claudia. I. Henschke, D. P. 
Naidich, D. F. Yankelevitz, G. McGuinness, D. I. McCauley, J. P. Smith, D. M, 
Libby, M. W. Pasmantier, M. Vazquez, J. Koizumi, D. FUeder, N. K. Altorki, and O. 
S. Miettinen, "Early Lung Cancer Action Project: Initial Findings on Repeat 
5 Screening." Cancer July 1 , 200 1 ; 92(1 ): 1 53- 1 59. The introduction of Computed 
Tomography (CT) scanners, particularly scanners with helical capabilities, has 
increased the resolution of lung images and greatly increased the number of images 
per screening study that must be evaluated by the radiologist. The development of 
the computed tomography (CT) technology and post-processing algorithms has 

10 provided radiologists with a useful tool for diagnosing lung cancers at early stages. 
However, current CT systems have their inherent shortcomings in that the amount of 
chest CT images (data) that is generated from a single CT examination, which can 
range from 30 to over 300 slices depending on image resolution along the scan axial 
direction, becomes a huge hurdle for the radiologists to interpret. Accordingly, there 

15 is a constant need for the improvement and development of diagnostic tools for 

enabling a radiologist to review and interpret the vast amount of information that is 
obtained through a CT examination. 

International Publication No. WO 01/78005 A2 discloses a system and method 
20 for three dimensional image rendering and analysis, and is incorporated herein by 

reference. The system performs a variety of tasks that aid a radiologist in interpreting 
the results of a CT examination. 

U.S. Patent Application Serial No. 10/245,782 discloses a system and method 
25 directed to diagnostic imaging of small pulmonary nodules, and is incorporated herein 
by reference. The application includes methods for detection and feature extraction 
for size characterization, and focuses on the analysis of small pulmonary nodules that 
are less than 1 centimeter in size, but is also suitable for larger nodules as well. 

30 Radiologist generally fail to detect nodules primarily due to interpretation 

and oversight error. The use of a computer aided detection (CAD) systems avoids 
these human related errors to tremendously improve diagnostic accuracy. M. Fiebich, 
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C. Wietholt, B.C. Renger, S.G Armato, K.R. Hoffimann, D. Wormanns and S. 
Diederich, "Automatic detection of pulmonary nodules in low-dose screening thoracic 
CT ex-aminations", SPIE, vol. 3661, pp 1434-1439, 1999; S.G Armato, 1 1 1, F. Li, 
M.L. Giger, "Performance of Automated CT Lung Nodule Detection on Missed 
5 Cancers", scientific paper presentation, RSNA 87th scientific assembly and aimual 
meeting, Nov 25 -30, 2001; F. Li, S. Sone, H. Abe, H.M MacMahon, S.G Armato, K. 
Doi, "Missed Lung Cancers in low-dose Helical CT Screening Obtained from a 
General Population", scientific paper presentation, RSNA 87th scientific assembly 
and annual meeting, Nov 25 -30, 2001 ; C.L Novak, D.P Naidich, L. Fan, J. Qian, J.P. 
10 Ko, A.N Rubinowitz, "Improving Radiologists* Confidence of Interpreting Low-dose 
Multidetector Lung CT screening Studies Using an Interactive CAD system". 
Scientific paper presentation, RSNA 87th scientific assembly and atmual meeting, 
Nov 25 -30, 2001 . CAD systems perform automated nodule detection in addition to 
providing useful visualization tools for the radiologists. 

15 

Novak et al reported an improvement in detection of potential nodules from 22 
to 77% with the availability of an interactive CAD (ICAD) system for the radiologist. 
C.L Novak, D.P Naidich, L. Fan, J. Qian, J.P. Ko, A.N Rubinowitz, "Improving 
Radiologists* Confidence of Interpreting Low-dose Multidetector Lung CT screening 

20 Studies Using an Interactive CAD system". Scientific paper presentation, RSNA 87th 
scientific assembly and annual meeting, Nov 25 -30, 2001 . Potential nodules were 
identified and rated with and without the ICAD tools for nodule interpretation. They 
concluded from their results that the interactive CAD systems significantly increase 
radiologists confidence when interpreting CT screening studies. In another study, 

25 Aramato et al reported that 78% of nodules missed during visual interpretation were 
detected by their automated method. 72% of the missed nodules were due to 
oversight* error, and the rest were due to interpretation error. S.G Armato,l 1 1, F. Li, 
M.L. Giger, "Performance of Automated CT Lung Nodule Detection on Missed 
Cancers", scientific paper presentation, RSNA 87th scientific assembly and annual 

30 meeting, Nov 25 -30, 2001 . All the nodules missed due to oversight error were 
detected by their CAD system whereas 70% of the nodules missed due to 
interpretation error were detected. Li et al also showed that 78% of missed nodules 



were detected by their computerized scheme. F. Li, S. Sone, H. Abe, H.M 
MacMahon, S.G Armato, K. Doi, "Missed Lung Cancers in low-dose Helical CT 
Screening Obtained from a General Population", scientific paper presentation, RSNA 
87th scientific assembly and aimual meeting, Nov 25 -30, 2001. According to them, 
5 the main reason for detection errors were difficulty in detection due to small size and 
low intensity, oversight due to adjacent or overlapping pulmonary vessels or fissures 
and lack of attention to relatively obvious nodules adjacent to the pulmonary hilum. 
In a related study, Fiebich et al reported a 15% improvement in detection sensitivity 
using their CAD system in addition the conventional film reading procedure, M. 
10 Fiebich, C. Wietholt, B.C. Renger, S.G Armato, K.R. Hoffinaim, D. Wormanns and S. 
Diederich, "Automatic detection of pulmonary nodules in low-dose screening thoracic 
CT ex-aminations", SPIE, vol, 3661, pp 1434-1439, 1999. 

The evolution of CT scaimer technology has played an important role in the 
15 development of detection algorithms. Early low resolution whole lung CT scans had 
a slice thickness of 5-1 0mm with a 0.5-0.6 mm in-plane resolution. Because of this 
low axial resolution, early computer detection algorithms were based entirely on two 
dimensional (slice-by-slice) image analysis techniques. Currently, multi-slice helical 
scanners with better axial resolution are widely available. This improvement in axial 
20 resolution permits three dimensional image analysis techniques which can detect 
smaller nodules. 

Recent research on pulmonary nodule detection has focused on 3D region 
identification and feature extraction procedures followed by rule-based classification. 

25 Aramato et al implemented a computerized scheme that uses 2D and 3D extracted 
features from regions identified by multiple level gray-level thresholding. S.G. 
Armato III, M. L. Giger, J.T Blackburn, K. Doi, H. MacMahon, "Three-dimensional 
approach to lung nodule detection in helical CT", SPIE, vol. 3661, pp. 553-559, 1999. 
In this paper, they used a rolling ball algorithm to avoid missing nodules attached to 

30 the pleural surface. They reported an operating point of 85% sensitivity and 89% 
specificity indicating an overall sensitivity of 70% with an average of three false- 
positive per slice. Similarly, 2D and 3D geometrical features have been used by 



Gurcan et al in their detection algorithm. M.K. Gurcan, N. Petrick, B. Sahiner, H.P. 
Chan, P.N. Cascade, E.A. Kazerooni, L.M. Hadjiiski, "Computerized lung nodule 
detection on thoracic CT images: combined rule-based and statistical classifier for 
false positive reduction", SPIE, Vol. 4322, pp 686-692, 2001 . They reported a 84% 
5 detection rate with 1.75 FPs per slice detection results tested on 17 patients with a 
total of 3 1 lung nodules. Fan et al implemented an adaptive 3D region growing 
algorithm followed by a classification scheme that makes use of geometric features 
such as diameter, volume, sphericity, mean intensity value and standard deviation of 
intensity. L. Fan, C. Novak, J. Qian, G. Kohl and D. Naidich, "Automatic Detection 

10 of Lung Nodules from Multi-Slice Low-Dose CT Images", SPIE, vol.4322, pp 1828- 
1835, 2001. This algorithm only detects nodules with very small vasculature 
connections and no large solid structure attachment. Toshioka et al tested their 
detection algorithm which on 450 cases (15,750 images). S. Toshioka, K. Kanazawa, 
N. Niki, H. Satoh, H. Ohmatsu, K. Eguchi, et al, "Computer aided diagnosis system 

1 5 for lung cancer based on helical CT images", SPIE, vol. 3034, pp 975-984, 1 997. 
Compared with image interpretation by 3 Radiologists, CAD detected all tumors 
identified as highly probable with 5 false negatives (4 of which represented tumors 
less than 5mm in size) and 1 1 false positive cases(ranging firom 2.4/case for "high 
probability" nodules to 7.2/case for "suspicious" nodules). 

20 

Lee et al used a template matching technique based on a genetic algorithm to 
detect nodules. Different templates were generated for nodules with and without an 
attachment to the pleural surface. Y. Lee, T. Hara, H. Fujita, S. Itoh and T. Ishigaki, 
"Automated Detection of Pulmonary Nodules in Helical CT images Based on an 

25 Improved Template-Matching Technique", IEEE Transactions on Medical Imaging, 
Vol. 20, No. 7, pp 595-604, 2001. Although, they developed an elegant mathematical 
model of a nodule, the algorithm resulted in a very high number of false positives (4.4 
per slice) with 72% sensitivity. Other computer vison methods have also been 
explored for pulmonary nodule detection. Morphological analysis techniques have 

30 been utilized for detection of suspicious regions. H. Taguchi, Y. Kawata and N. Niki, 
H. Satoh, H. Ohmatsu, K. Eguchi, M. Kaneko and N. Moriyama, "Lung cancer 
detection based on helical CT images using curved surface morphology analysis". 
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SPIE, vol 3661, pp 1307-1313, 1999. Penedo et al have developed a computer aided 
detection system based on a two level artificial neural network. M.G. Penedo, M.J. 
Carreira, A. Mosquera, and D. Cabello, "Computer-aided diagnosis: a neuralnetwork 
based approach to lung nodule detection", IEEE Transactions on Medical Imaging, 
5 vol.17, no.6, pp 872-879, 1998. The first network performs detection of suspicious 
regions, while the second one classify the regions based on the curvature peak on all 
points in the suspicious region. They recorded 89% - 96% sensitivity with 5-7 FPs 
per slice. Artificial neural networks* capabilities have also been used by Lo et el. S-C 
B. Lo, S-L A. Lou, J-S Lin, M.T. Freedman, M V. Chien and S. K, Mun, "Artifical 
10 convolution neural Network techniques and applications for lung nodule detection", - 
IEEE transactions on medical imaging, vol. 14, no.4, pp 711-718, 1995. In their 
work, Lo et al used a convolution type neural network which recorded an 82% 
detection rate. 

1 5 Object-based deformation techniques have been incorporated into detection 

systems. Lou et al used deformation techniques to differentiate lung nodules from 
blood vessels in their 3D CT lung nodule detection system. S.L Lou, C.L Chang, 
K.P Lin and T. Chen, "Object based deformation technique for 3-D CT lung nodule 
detection", SPIE, vol 3661, pp 1544-1552.1999. This research did not address 

20 surface irregularities that occur in nodules with significant vasculature connections. 
Knowledge based techniques have also been used in recent research. Works of 
Erberich et al and Brown et al can be cited in this regard. S.G. Erberich, K.S. Song, 
H. Arakawa, H.K Huang, R. Webb, K. S. Hoo, B.W. Loo, "knowledgebased lung 
nodule detection from helical CT", RSNA 1997 annual meeting; M.S. Brown and M. 

25 F. McNitt-Gray, "Method for segmenting chest CT image data using an anatomical 
model: preliminary results", SPIE, vol- 16, no.6, pp 828-839, 1997. Erberich et al 
used rule based tree to classify candidates generated using gradient Hough 
transformation. Detection performance statistical results were not reported in their 
paper. Brown et al developed a multipurpose modular knowledge based system. 

30 They demonstrated nodule detection application using this modular architecture. 
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A review of the prior art indicates that progress has been made on 
developing automated detection programs for lung nodules in helical CT scans. 
However, there is a large variations in performance, likely caused by the small data 
sets used in these studies. Much more effort is need to bring the performance of these 
5 computerized detection systems to level acceptable for clinical implementation. Most 
of the detection algorithms have been designed to detect a single type of a nodule (i.e 
nodule with a small vessel coimection). Nodules with significant vessel coimections 
or attachment to large solid structure have been either reported as a missed or not 
considered in the detection performance evaluation. Y. Lee, T. Hara, H, Fujita, S. Itoh 

10 and T. Ishigaki, "Automated Detection of Pulmonary Nodules in Helical CT images 
Based on an Improved Template-Matching Technique", IEEE Transactions on 
Medical Imaging, Vol.20, No.7, pp 595-604, 2001; L. Fan, C. Novak, J. Qian, G. 
Kohl and D. Naidich, "Automatic Detection of Lung Nodules from Multi-Slice Low- 
Dose CT Images", SPIE, vol.4322, pp 1828-1835, 2001; M. Fiebich, C. Wietholt, 

15 B.C. Renger, S.G Armato, K.R. Hoffmann, D. Wormanns and S. Diederich, 

"Automatic detection of pulmonary nodules in low-dose screening thoracic CT ex- 
aminations", SPIE, vol. 3661, pp 1434-1439, 1999. Accordingly, there is a great need 
for an algorithm which detects nodules with or without attachments to large solid 
structures with fewer false positives. 

20 

One predictor of malignancy of a pulmonary nodule in a CT image is the 
change in volume of the nodule. The change in volume can be measured as percent 
volume change or a doubling time. To obtain these measurements, two high- 
resolution CT scans, separated by a few months, are taken of the nodule. The 

25 nodules are segmented from the CT images and the percent volume change or 

doubling time is calculated using the segmented nodule volumes. The accuracy of 
the change in volume measurement is dependent on the consistency of the 
segmentations of the nodule in the two images. In the extreme case, a 
missegmentation of one of the nodules may adversely affect the malignancy 

30 predictor by moving the doubling time measurement above or below the threshold 
for malignancy. 
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There has been some work on tracking the change of pulmonary nodules in 
CT images. In Kawata et al, the pulmonary nodules are registered together using 
rigid-body registration and affine registration at two different stages. Y. Kawata, N. 
Niki, H. Ohmatsu, M. Kusumoto, R. Kakinuma, K. Mori, N. Nishiyama, K. Eguchi, 
5 M. Kaneko, and N. Moriyama. Tracking interval changes of pulmonary nodules using 
a sequence of three-dimensional thoracic images. In Medical Imaging 2000: Image 
Processing, Proceedings ofSPIE, volume 3979, pages 86-96, 2000. The nodules are 
segmented using a 3-D deformable surface model and curvature features are 
calculated to track the temporal evolution of the nodule. This work was extended by 

10 Kawata et al, by adding an additional 3-D non-rigid deformable registration stage and 
the analysis was performed using a displacement field to quantify tiie areas of nodule 
growth over time. Y. Kawata, N. Niki, H. Ohmatsu, M. Kusumoto, R. Kakinuma, K. 
Mori, N. Nishiyama, K. Eguchi, M. Kaneko, and N. Moriyama. Analysis of evolving 
processes in pulmonary nodules using a sequence of three-dimensional thoracic 

15 images. In M. Sonka and K. M. Hanson, editors. Medical Imaging 2001: Image 

Processing, Proceedings ofSPIE, volume 4322, pages 1 890- 1 90 1 , 200 1 . In Reeves 
et al, a method was introduced to estimate the growth of a nodule without the explicit 
use of segmentation. A. P. Reeves, W. J. Kostis, D. F. Yankelevitz, and C. I. 
Henschke. Analysis of small pulmonary nodules without explicit segmentation of CT 

20 imsigQS. Radiology, 21 7P:243-244, November 2000. The pulmonary nodules are 
registered using translation and the doubling time is calculated from the gaussian- 
weighted regions-of-interest. In Kostis et al, and Reeves et al, a segmentation method 
based on thresholding and morphological filtering is discussed. W. J. Kostis, A. P. 
Reeves, D. F. Yankelevitz, and C. I. Henschke. Three-dimensional segmentation of 

25 solitary pulmonary nodules from helical CT scans. In H. U. Lemke, M. W. Vannier, 
K. Inamura, and A. G. Farman, editors. Proceedings of Computer Assisted 
Radiology and Surgery (CARS '99), pages 203-207. Elsevier Science, June 1999; A, 
P. Reeves and W. J. Kostis. Computer-aided diagnosis of small pulmonary nodules. 
Seminars in Ultrasound, CT, and MRI, 21(2): 116-128, April 2000. From the 

30 nodule segmentation, the volume of the nodule can be easily calculated and the 
doubling-time can be determined. 



SUMMARY OF THE EVVENTION 

The present invention is directed to diagnostic imaging of small pulmonary 
5 nodules. There are two main stages in the evaluation of pulmonary nodules from CT 
scans: detection, in which the locations of possible nodules are identified, and 
characterization, in which a nodule is represented by measured features that may be 
used to evaluate the probability that the nodule is cancer. Currently, the most usefril 
prediction feature is growth rate, which requires the comparison of size estimates 
10 from two CT scans recorded at different times. The present invention includes 

methods for detection and feature extraction for size characterization. The invention 
focuses on small pulmonary nodules that are less than 1 centimeter in size, but these 
methods are also suitable for larger nodules as well. 

1 5 For the purpose of Computer Aided Diagnosis (CAD), Pulmonary nodules 

are dichotomized into attached nodules and isolated nodules based on their 
location with respect to other solid lung structures. Attached nodules are adjacent 
to some larger solid structure, such as the pleural surface. Isolated nodules consist 
of both well-circumscribed nodules and nodules that are larger than all adjacent 

20 structures, such as blood vessels or bronchi. The nodules may be solid, non-solid 
or part-solid. The overall analysis of a Computed Tomography scan generally 
entails the following: 

1. Detection 

25 (a) Identify the lung regions and main bronchi from thoracic CT images 

(b) Separate the lungs into two major regions: (1) the lung parenchyma 
and (2) the lung surface region, including the pleural surface and major airways. 

(c) Identify possible locations of isolated nodules in the lung parenchyma 
region and identify possible locations of attached nodules in the in the lung surface 

30 regions. 
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2. Characterization 

(a) Starting with a single location point within a possible nodule, 
identify the nodule region in the CT images. This entails locating the geometric 
center of the nodule and approximating its size. 
5 (b) Given the location and approximate size of a nodule, compute 

characteristic features of the nodule, including robust size estimates. 

The present invention identifies possible locations of isolated nodules in the 
lung parenchyma region within whole lung CT scans. With reference to the overall 

10 analysis listed above, this aspect corresponds to the first part of task 1(c). In 

particular, the invention addresses the detection problem of isolated nodules and 
expands upon the lung segmentation algorithms disclosed in United States Patent 
application serial no. 10/245,782. The invention searches the segmented lung region 
with a filter to identify lung nodule candidates. A set of filters then eliminate a 

1 5 substantial amount of false candidates until only those with a high probability of 
being true nodules remain. 

The present invention also includes methods for improving the consistency of 
pulmonary nodule segmentation using 3D rigid-body registration. The accurate 
20 repeat segmentation of lung nodules is critical to the measurement of nodule growth 
rate, which is the most accurate non-invasive computer based predictor of 
malignancy. The invention increases the accuracy of repeat segmentation by image 
comparison methods including (a) 3D rigid body registration, (b) histogram matching, 
and (c) knowledge-based matching of the vessel removal. 

25 

A preferred embodiment of the present invention is a method and apparatus 
for analyzing a computed tomography scan of a whole lung for lung nodules. The 
apparatus of the invention is a detecting unit configured with the teachings of the 
method of the invention. The invention can also be practiced on a machine readable 
30 medium. The method includes segmenting a first lung region and a second lung 
region from the computed tomography scan. The first lung region corresponds to 
lung parenchyma of the lung and the second lung region corresponds to at least one of 

10 
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a pleural surface of the lung and a surface defined by vessels within the lung. An 
initial list of nodule candidates is generated from the computed tomography scan 
within the first lung region. The list includes at least a center location and an 
estimated size associated with each nodule candidate. Subimages are next generated 
for each nodule candidate in the initial list. Streaking artifacts are then selectively 
removed from the subimages. The nodule candidates are filtered to eliminate false 
positives from the list. 

In a preferred embodiment, the initial list is generated by thresholding the first 
lung region. Nodule candidate regions are next identified by labeling high density 



voxels foreground voxels. R mi is determined for each foreground voxel. The local 



maximum R mi is selected within a nodule candidate region. The limited extent 



criterion is next determined for each foreground voxel which corresponds to a Rmi - 
The initial list of nodule candidates is generated for nodules which satisfy the limited 



extent criterion. The list includes at least Nc and R mi associated with the 
corresponding foreground voxel. 

In another preferred embodiment, the streaking artifact removal is 
accomplished by initially determining the amount of streaking artifact present in the 
sub-image followed by filtering the streaking artifact out from the subimage when the 
amount of the streaking artifact present in the sub-image exceeds Tsar. Preferably the 
amount of streaking artifact present in the sub-image is calculated by a metric 



Preferably the filtering is performed by a vertical median filter of size 1x3 and Tsar is 
selected to be in a range from about 20000 to about 80000. 
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In another preferred embodiment, the false positives are eliminated from the 
list by initially determining for each nodule candidate a fraction. Fay of a surface of 

11 



the nodule candidate that is attached to other solid structures followed by removing 
the nodule candidate from the list when the fraction exceeds Ta. 

In another preferred embodiment, the false positives are eliminated from the 
list by initially generating a cube wall about each nodule candidate followed by 
determining an intersection volume, Vn^ corresponding to portions of the nodule 
region associated with the nodule candidate that intersect the cube wall. Nodule 
candidates are removed from the list when the fraction of the intersection volume, Vnu 
over the volume of the nodule candidate, exceeds 7Vv. 

A preferred embodiment of the present invention is a method and apparatus 
for correlating a segmentation of 3-d images of a pulmonary nodule from a high- 
resolution computed tomography (CT) scans. The images include a first image {im\) 
obtained at time-1 and a second image {imi) obtained at time-2. The apparatus of the 
invention is a registration unit configured with the teachings of the method of the 
invention. The invention can also be practiced on a machine readable medium. The 
method includes selecting a first region-of-interest (ROIi) for the nodule in the first 
image (zmi) and selecting a second region-of-interest (ROI2) for the nodule in the 
second image (//W2). The second region-of-interest (ROI2) is registered to the first 
region-of-interest (ROIi) to obtain a transformed second region-of-interest (ROIzr). 
The nodule in the first region-of-interest (ROIi) and the transformed second region- 
of-interest (ROI2O are both separately segmented. The first segmented nodule 
and the second segmented nodule (52) are then adjusted. 

In a preferred embodiment, the nodule in the first region-of-interest (ROIi) 
and the transformed second region-of-interest (ROtr) are both separately segmented 
by performing at least one of the following: 

(i) gray-level thresholding; 

(ii) morphological filtering for vessel removal: and 

(iii) plane clipping for separarting a pleural wall. 
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Preferably the gray-level thresholding is performed at an adaptive threshold level. 
The adaptive threshold level is preferably selected for each region-of-interest (ROIi 
and ROI2) by: 

determining a peak parenchyma value, Vp\ 
determining a peak nodule value, v„; 

calculating the adaptive threshold level as a midpoint between the peak 
parenchyma value, Vp, and the peak nodule value, v^. 
An intensity histogram, H{pC) is preferably calculated for determining the peak 
parenchyma value, Vp, and the peak nodule value, v„. The intensity histogram, H{x)^ is 
preferably filtered with a gaussian with standard deviation of about 25 HU prior to 
determining the peak parenchyma value, Vp, and the peak nodule value, v^. 

In another preferred embodiment, the registration to obtain a transformed 
second region-of-interest (ROI2O is performed by: 

(a) calculating initial rigid-body transformation parameters for a 
rigid-body transformation on the second region-of-interest (ROI2); 

(b) determining the optimum rigid-body transformation parameters 
by calculating a registration metric between the first region-of-interest (ROIi) 
and the rigid-body transformation on the second region-of-interest (ROI2); and 

(c) generating a registered image from the optimum rigid-body 
transformation parameters. 

Preferably the registration metric is calculated by: 

transforming the second region-of-interest (ROI2) witti the initial rigid- 
body transformation parameters to obtain a transformed second region-of- 
interest (ROl2f); 

calculating the registration metric as a mean-squared-difference (MSD) 
between the transformed second region-of-interest (ROI2O and the first region- 
of-interest (ROIi); and 

searching for the minimum mean-squared-difference (MSD) in the 6- 
dimensional parameter space. 
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The transforniing of the second region-of-interest (ROI2) to obtain the transformed 
second region-of-interest (ROI2/) is preferably a mapping of a point v in 3-d space to a 
point v' in transformed space defined by : 

tx 



V=RxRyRzv^ 



ty 

tz 



wherein Rx^ Ry, and Rz are rotation matrices defined as: 

10 0 



R. = 



R. = 



R. = 



0 cos(r^) -sin(r^) 
0 sin(rj cos(rJ 

cos(r^) 0 sin(r^) 

0 1 0 
-sin(r ) 0 cos(r ) 



cos(r, ) 
sin(rj 
0 



-sin(rj 0 
cos{r^ ) 0 
0 1 



10 The initial rigid-body transformation parameters preferably include six parameters 

(tx^ty^tz^rx^ry^rz) respectively defined as translation in x, translation in y, translation in 
z, rotation about the x-axis, rotation about the y-axis, and rotation about the z-axis. 
The initial rotation parameters (rx^ry^rz) are generally all set to zero while the initial 
translation parameters (pc,ty,tz,) are set so that the nodule in the first region-of-interest 

1 5 (ROIi) overlaps the nodule in the second region-of-interest (ROI2) during the initial 
calculation of the registration metric. Preferably the mean-squared-difference (MSD) 
is gaussian weighted. 

In another preferred embodiment, a first thresholded image (ri)and a second 
20 thresholded image (T2) are defined by gray-level thresholding prior to vessel removal 
and separating the pleural wall. The adjustment of first segmented nodule (S\) and the 
second segmented nodule (52) is performed by comparing the segmented nodules and 
the thresholded images. The active pixels in the segmented nodules are marked as 
one of: 

25 a repeat nodule pixel; 
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a nodule growth pixel; 

a nodule atrophy pixel; and 

a nodule missegmentation pixel. 



5 A foreground pixel in the first segmented nodule is preferably marked as a 
repeated nodule pixel from the first region-of-interest (ROIi) to the transformed 
second region-of-interest (ROI2/) when the corresponding pixel in second segmented 
nodule (^2) and the corresponding pixel in second thresholded image (Tz) are both 
foreground. A foreground pixel in the first segmented nodule (S\) is preferably 

10 marked as a nodule atrophy pixel when the corresponding pixel in second segmented 
nodule (Si) is background and the corresponding pixel in second thresholded image 
(T2) is background. A foreground pixel in the first segmented nodule (S\) is 
preferably marked as a missegmented pixel in the first region-of-interest (ROIi) when 
the corresponding pixel in second segmented nodule (52) is background and the 

15 corresponding pixel in second thresholded image (T2) is foreground. A foreground 
pixel in the second segmented nodule (S2) is preferably marked as a repeated nodule 
pixel from the first region-of-interest (ROIi) to the transformed second region-of- 
interest (ROI2/) when the corresponding pixel in first segmented nodule (Si) and the 
corresponding pixel in first thresholded image (T\) are both foreground. A foreground 

20 pixel in the second segmented nodule (52) is preferably marked as a nodule growth 
pixel when the corresponding pixel in first segmented nodule (5i) is background and 
the corresponding pixel in first thresholded image {T\) is background. A foreground 
pixel in the second segmented nodule (^2) is preferably marked as a missegmented 
pixel in the transformed second region-of-interest (ROI2O when the corresponding 

25 pixel in first segmented nodule (5i) is background and the corresponding pixel in first 
thresholded image (T\) is foreground. 

For a better understanding of the present invention, reference is made to the 
following description to be taken in conjunction with the accompanying drawings and 
30 its scope will be pointed out in the appended claims. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

Preferred embodiments of the invention have been chosen for purposes of 
illustration and description and are shown in the accompanying drawings, wherein: 

5 

Figure 1 is a flowchart of the automated detection algorithm; 
Figure 2 is an original single axial image scan; 

Figure 3 illustrates the lung parenchyma after being segmented from the image 
shown in Figure 2; 

10 Figure 4 illustrates the pleural surface corresponding to the surface of the 

volume shown in Figure 3; 

Figure 5 illustrates an isolated pulmonary nodule; 

Figure 6 illustrates an pulmonary nodule attached to pleural surface; 

Figure 7 is a 2-dimensional illustration of Rm and Rmc \ 
1 5 Figure 8 is a 2-dimensional illustration of volume occupancy procedure; 

Figure 9 illustrates the spherical shell used for nodule adjacent region analysis; 

Figure 10 illustrates a nodule sub-image before streaking artifact removal; 

Figure 1 1 illustrates a nodule sub-image after streaking artifact removal; 

Figure 12 illustrates a multi-stage filtering procedure for successive 
20 refinement of the nodule candidates; 

Figure 13 illustrates nodule candidate attachment analysis; 

Figure 14 illustrates the region growing procedure for the attachment analysis; 

Figure 15 illustrates the effect of global thresholding on small vessel 
bifiircation points on the original gray scale image; 
25 Figure 16 illustrates the effect of global thresholding on small vessel 

bifiircation points on the binary image; 

Figure 17 illustrates small bifiircation point removal; 

Figure 18 is a histogram of a real 12mm nodule between -1200HU and 200 

HU; 

30 Figure 19 is a histogram of a real 12mm nodule between -200HU and 200HU; 

Figure 20 is a histogram of an ideal nodule between -1200HU and 200 HU; 
Figure 21 is a histogram of an ideal nodule between -1200HU and 200 HU; 
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Figure 22 is a histogram of a real 5mm nodule between -1200HU and 200HU; 
Figure 23 is a histogram of a real 5mm nodule between -200HU and 200HU; 
Figure 24 illustrates an image of the rigid-body registration of the first ROI; 
Figure 25 illustrates an image of the rigid-body registration of the second ROI; 
5 Figure 26 illustrates an image of the rigid-body registration of the second ROI 

registered to the first ROI; 

Figure 27 illustrates the difference image between the first ROI and the second 

ROI; 

Figure 28 illustrates the difference image between the first ROI and the 
1 0 registered second ROI; 

Figure 29 illustrates a montage of the original gray scale image showing the 
nodule region-of-interest; 

Figure 30 illustrates a montage showing the segmented nodule corresponding 
to Figure 31; 

1 5 Figure 3 1 illustrates a first thresholded image (T\) corresponding to the first 

region-of-interest (ROIj); 

Figure 32 illustrates a second thresholded image (72) corresponding to the 
transformed second region-of-interest (ROl2^); 

Figure 33 illustrates a first segmented nodule (Si) corresponding to the first 
20 region-of-interest (ROI j ); 

Figure 34 illustrates a second segmented nodule (52) corresponding to the 
transformed second region-of-interest (ROl2^); 

Figure 35 illustrates the first segmented nodule (Si) having the active pixels 
therein marked as (A) nodule in time 1 , (B) nodule atrophy, or (C) missegmentation; 
25 Figure 36 illustrates a second segmented nodule (^2) having the active pixels 

therein marked as (D) nodule in time 2, (E) nodule growth, or (F) missegmentation; 

Figure 37 illustrates a first adjusted segmented nodule (N\) corresponding to 
the first region-of-interest (ROIj); and 

Figure 38 illustrates a second adjusted segmented nodule (A^2) corresponding 
30 to the transformed second region-of-interest (ROI2/). 
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DETAILED DESCRIPTTON OF THF. INVENTION 



A system in accordance with the present invention may include a scanner, 
processor, memory, display device, input devices, such as a mouse and keyboard, 
and a bus connecting the various components together. The system may be coupled 
to a communication medium, such as a modem connected to a phone line, wireless 
network, or the Internet. 

The present invention is preferably implemented using a general purpose 
digital computer, microprocessor, microcontroller, or digital signal processor 
programmed in accordance with the teachings of the present specification, as will be 
apparent to those skilled in the computer art. Appropriate software coding may be 
readily be prepared by skilled programmers based on the teachings of the present 
disclosure, as will be apparent to those skilled in the software art. 

The present invention preferably includes a computer program product, which 
includes a storage medium comprising instructions that can be used to direct a 
computer to perform processes in accordance with the invention. The storage 
medium preferably includes, but is not limited to, any type of disk including floppy 
disks, optical data carriers, compact discs (CD), digital video discs (DVD), magneto- 
optical disks, read only memory (ROM), random access memory (RAM), electically 
programmable read only memory (EPROM), electrically eraseable programmable 
read only memory (EEPROM), magnetic or optical cards, or any type of media 
suitable for storing information. 

Stored on any one of the above described storage media, the present invention 
preferably includes programming for controlling both the hardware of the computer 
and enabling the computer to interact with a human user. Such programming may 
include, but is not limited to, software for implementation of device drivers, operating 
systems, and user applications. Such storage media preferably fiirther includes 
programming or software instructions to direct the general purpose computer to 
perform tasks in accordance with the present invention. 
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' The programming of the computer preferably includes software for digitizing 
and storing images obtained from the image acquisition device (helical computed 
tomography scanner). Alternatively, it should be understood that the present 
invention may also be implemented to process digital data derived from images 
obtained by other means, such as x-rays and magnetic resonance imaging (MRI), 
positron emission tomography (PET), ultrasound, optical tomography, and electrical 
impedance tomography. 

The invention may also be implemented by the preparation of application 
specific integrated circuits (ASIC), field programmable gate arrays (FPGA), or by 
interconnecting the appropriate component devices, circuits, or modules, as will be 
apparent to those skilled in the art. 

A. NODULE DETECTION 

Referring to Figure 1, the nodule detection algorithm involves four major 
stages. First, lung regions are segmented from the whole lung computed 
tomography (CT) scan. This is followed by a hypothesis generation stage. In this 
stage, nodule candidate regions are identified from the whole lung scan and their 
size is estimated. In the third stage, the nodule candidate sub-images pass through a 
streaking artifact removal process. The nodule candidates are successively refined 
in the fourth stage using filters of increasing complexity. 

The region of the image space where pulmonary nodules are to be found is 
first identified. A distinction is made between the lung parenchyma where spherical 
nodules are located and the region of lung parenchyma adjacent to solid structures 
where attached nodules are located, since different techniques are used to detect the 
two nodule forms. Therefore, two lung regions are identified; the first lung region is 
the lung parenchyma which is not close to any major solid structure, and the second 
lung region is a narrow region of lung parenchyma that is adjacent to solid structures. 
In this context, major solid structures include the chest wall, the hilum region, large 
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blood vessels and the primary bronchi. An axial slice image and its corresponding 
segmented regions are shown in Figures 2 through 4. 

Pulmonary solid nodules are approximately spherical shaped lesions with a 
5 density slightly greater than water. These lesions are identified on routine chest 
radiographs and/or low resolution CT scans of asymptomatic patients. Pulmonary 
nodules can be broadly classified into two groups. The first group include nodules 
with no attachment to a large solid structure These nodules have, in general, a 
spherical shape. However, the spherical shape may be distorted by other small lung 
10 structures such as vessels, bronchi, scars and regions of morbidity. Typically, 
growing nodules appear to wrap around such structures. An example of such a 
nodule is shown in Figure 5. 

The second group consists of nodules attached to large dense structures, an 
1 5 example of which is shown in Figure 6. In this context, large dense structures refer to 
those structures with a size comparable to or larger than the nodules diameter. 
Nodules attached to the pleural surface are the most prevalent of this type. However, 
nodules may also be attached to other structures such as large vessels, airways and/or 
the hilum region. In general, the shape of this type of nodule is significantly affected 
20 by a large adjacent structure; a nodule that is attached to a pleural surface often has 
the appearance of a hemisphere with the curved region growing into the lung 
parenchyma. The images in Figures 5 and 6 were created using light shaded 
rendering of the surface of the nodule extracted from a high resolution CT image 
series. 

25 

The current radiological size characterization of a nodule is its "diameter" Nr. 
For CT image scans, this diameter is estimated from a single axial slice from the 
central region of the nodule. One common method of size estimation is computing 
the average of the horizontal and vertical extents of the nodule. 

30 

For three-dimensional computer image analysis, two parameters are defined: 
the maximum inscribed sphere (MIS) and the minimum circumscribed sphere (MCS) 

20 
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that characterize the size and irregularity of a nodule in three-dimensional space. The 
MIS and MCS are defined with respect to the geometric center of mass (GCOM). 
The GCOM Nc = (Xc Yc, Z^) is the location of the center of mass of the nodule region 
assuming a uniform density (I). The GCOM for a nodule is defined as follows: 

Let the three dimensional region of the nodule be defined by fNix,y,z). 



1 , if(x,y,z) lies within the nodule region 
^ • otherwise 



Then, the coordinates of the GCOM Nc=(Xc, Yc, Zc) are given by: 



^ jjjxfN(x,y,z)dxdydz 

J J J M^^y^^) ^ 



(1) 



Yc = ^f. (2) 

J J J fN{x,y,z)dxdydz 

I I I z fN{x, y, z) dx dy dz 

Zc = ^, (3) 

J J J Mx,y,z)dxdydz 



This measurement does not consider density variations within the nodule but 
depends only on its geometric form. For a nodule, the minimum distance (i?M/;from 
the nodule geometric center (Nc,) to the boundary of the nodule (Nb), specifies the 
25 radius of the MIS. That is, the maximum sized sphere that can be inscribed within 
the nodule centered on 



_ min{D{Nc,N[{x,y,z))) 
Rmi - 
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where Nl denote the i'th boundary point of the nodule. D denotes the Euclidean 
distance given by 
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DEi[XuYMl^2j2.Z2]) = ^J{X^ - ^2)^ + (7. + Ylf + (Zl - Zlf (4) 

Similarly the maximum distance (Rmc) from Nc to the boundary of the 
nodule (Nb) specifies the radius of the MCS; i.e. the minimum sized sphere that 
can circumscribe the nodule centered on Nc touching the boundary of the nodule. 

_ max(D(Nc,NUx,y,z))) 
Rmc - y . 



For the case when a nodule has a highly irregular shape such that Nc does not lie 
within the nodule, the definition of Rmi is further refined to be the minimum distance 
to the boundary from Nc in which a change from inside to outside the nodule occurs. 
Similarly, the definition of Rmc can be modified to the maximum distance to the 
15 boundary from Nc in which a change from inside to outside the nodule occurs. 
However, such a situation will occur very rarely in practice. 

The concepts of Rmi and i?A/c are illustrated on a 2-dimensional example 
shown in Figure 7. In the Figure, Rmi and Rmc are the maximum inscribed and 
20 minimum circumscribed circle radii for the image region shown respectively. The 
extension to 3-dimension is straightforward. 

In the hypothesis generation stage, nodule candidates are identified from the 
whole lung scan and their size is estimated. It is critical that no nodules are missed 
25 at this stage. However, since the search space is very large, it is also important that 
this process be computationally efficient. The nodule candidates are refined in the 
subsequent stages. 

First, the lung parenchyma intensity is thresholded using a global threshold 
30 parameter (J^). This parameter can be determined empirically fi-om a training 

dataset to separate lung parenchyma tissue and high density solid structures. Tg is 
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preferably determined as approximately the mid-value between the mean value for 
solid tissue and the mean value for lung parenchyma. A preferred value for Tg is 
-574 HU. During the thresholding procedure, high density anatomical structure 
voxels having values higher than Tg are labeled as foreground voxels, creating a 
5 binary image. 

Let S(r,c) denote a spherical region with radius r centered at point c = 
(Xc, Zc). Consider a nodule N and its characteristic spheres (S(Rmi.Nc,) and 
S{Rmc9 Nc)) as illustrated in Figure 7. In general, for a nodule N, the center is 
10 always at Nc which may be dropped such that S(Rmi) =S(Rmi,Nc)^ A binary 
image representation of this nodule should have the following properties. 

a. All voxels within S(Rmi) should be one. 

b. All voxels outside S(Rmc) should be zero. 

15 c. Within the spherical shell region (S(Rmc) - S(Rmi))s there will be 

voxels of both values due to the irregular surface of the nodule. 

In order to generate a hypothesis for a nodule of size J? a// at a location P = (x, y, z) in 
the binary lung parenchyma image, the following two criteria be met: 
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a. "Solid center": A large fraction (Jy) of voxels in S(Rmi) centered at P 
must be one. 

b. "Limited extent": A large fraction (Tav) of voxels in the region 
(S(Rmc + 5) - S(Rmc)) centered at P must be zero. 

where 6 is a certain constant distance. The "Solid center" criterion(a) ensures that the 
nodule consists of a dense mass; the threshold Ty allows for some voxels to be zero 
due to image noise. 



30 While the nodule is bounded by radius Rm/, there may be other dense 

structures in the lung region some of which may be adjacent or touching the nodule 
N. The "Limited extent" criterion (b) verifies the finite extent of the nodule region. 
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However, it also allows a small amount of foreground voxels in the region 
immediately surrounding the nodule due to adjacent or attached structures. The 
preferred ranges and values for Tv, Tav and 6 as noted below were determined 
empirically from an analysis of a set of images in the training dataset. 

The hypothesis generation stage is implemented as follows: 



1. For each foreground voxel compute the maximum inscribed sphere 
radius, Rmi, for which the solid center criterion (a) is true. 
10 2. Select only the voxels with local maximum i?A47 values. 

3. For all voxels identified in step 2, determine the limited extent criterion 
(b). 

4. Make a list of nodule candidates of all voxels which satisfy the 
limited extent criterion (b) in step 3, including their center 

1 5 coordinates (iVc) and their size estimate Rmi. 



First, approximate spherical regions are identified and their size is 
estimated. Kanazawa et al. used gray weighted distance transforms for this 
purpose. K.Kanazawa, M.Kubo, N.Niki, H.Satoh, H.Ohmatsu, K.Eguchi, 

20 N.Moriyama, "Computer Assisted Diagnosis of Lung Cancer Using Helical X- 
ray CT", Proceedings of ICPR, pp 381-385, 1996. After computing the distance 
transform, Kanazawa et al. used a thinning procedure to determine the center point 
of the nodule candidate. In the present invention, a different noise tolerant 
procedure has been implemented. This procedure makes use of a volume 

25 occupancy calculation. Volume occupancy, at a particular foreground voxel (Fp), 
is defined as the ratio of the number of foreground voxels within a sphere {Sir^ 
Fp)) to the volume of the sphere. The procedure is illustrated for a 2D region in 
Figure 8 where Fp represents any seed pixel in the region R. For a 2D region, a 
circle is used rather than a sphere and an area coverage is calculated instead of a 

30 volume occupancy. The basic idea is to search for a center point pixel location in 
the region R at which the largest size circle can fit. This pixel corresponds to Nc 
for region R in Figure 8. 



The algorithm starts with an initial point Fp within the region and a circle of 
certain initial small radius. If the area coverage for this particular circle centered at 
Fp is greater than the threshold Jy the circle radius is incremented. This is 
5 sequentially done until the area coverage is less than Ty. Ty is preferably selected to 
be in a range between about .70 and 1.00, and most preferably about .80. The final 
circle size approximates the maximum circle size that can be fit in the region centered 
at location Fp. This procedure is performed for all the points in the region. Every 
point in the region would have a maximum circle radius associated with it. The 
10 region pixel that registered the largest circle size is marked as the center point and the 

A 

size is recorded as an estimate of the maximum inscribed circle ( R mi)- The 
computational efficiency of this procedure is significantly improved by updating the 
region pixels list at every iteration. The extension of this concept to 3D regions is 
straightforward using a spherical kemel and volume occupancy calculation as 
1 5 described above. 

After identifying the center point(A^c) of the nodule candidate region and 
measuring its maximum inscribed sphere's radius (Rmi), the next step is the 
determination and analysis of the immediately adjacent regions. 

20 

The volume occupancy of a spherical shell enclosing the nodule candidate is 
used to evaluate the degree of attached structures. The inner radius of this shell is 

A 

R MI i.e the estimated maximum inscribed sphere corresponding to Rmi- The middle 

A 

sphere has a radius of R mc which is obtained from a precomputed table of values for 
25 a given R^j^ The outer spherical shell has a preset thickness of 5 as shown in Figure 

A 

9. The R MC can also be calculated as the radius of a sphere having a volume between 
2.5 and 3.5 times that of the maximum inscribed sphere, Rmi- The preferred value of 

R MC can be calculated as \Fi*RMi. The 6 can be calculated as the radius of a sphere 
with a volume between 5.5 and 6.5 times the maximum inscribed sphere, Rml minus 

A ^ A 

30 R MC' The preferred value of 5 can be calculated as v6 *Rmi - R mc. 
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High density voxels in Region C correspond to attached structures. The center 
location Nc is then labeled as a nodule candidate if the ratio of the total volume (Vc) of 
high density voxels in the adjacent region C to the nodule candidate volume (Vn) is 
5 less than a certain threshold (Tav)- This rule is referred to as the adjacent region rule 
(AJR): 

Vc 

— < r.v (5) 

where Tav is the threshold value. Tav is preferably selected to be in a range between 
10 about 0.00 and 0.50, and most preferably about .28. 

After the initial list of nodule candidates is generated, a subimage for each 
candidate is preferably generated by clipping a region surrounding each nodule 
candidate to simplify data management. The subimages are basically segmented 
1 5 from the original whole lung CT scan. In an alternative embodiment of the 
invention, the subimage generation step could be skipped and the refinements 
described below (e.g. streaking artifact removal and multi-stage filtering) could be 
applied to the whole lung CT scan. 

20 The occurrence of CT image artifacts such as streaking poses a major 

problem for the detection of very small nodules. Streaking artifacts are caused by 
starvation of the x-ray photon flux and beam hardening effects. A majority of the 
streaking artifact occur near the patient shoulder area or when the patient arms are 
inside the scan FOV. Photon deficiency is limited to the projections that pass 

25 through both shoulders of the patient and result in a horizontal streaking pattern. 
Jiang Hsieh, "Generalized adaptive median filters and their application in 
computed tomography", SPIE, vol. 2298, pp. 662-669, 1994. 

Geometric characteristics of nodule candidates are used in the nodule 
30 candidate refinement stage. Artifacts deform geometrical properties of nodules 
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resulting in trae nodule elimination. A streaking artifact filter is preferably selectively 
applied to avoid deformation or elimination of small nodules. 

The present invention implements an adaptive streaking artifact removal 
5 filtering technique. After nodule candidate generation, individual nodule candidate 
sub-images are selectively filtered based on the amount of streaking artifact present. 
For this purpose, streaking artifact quantification metric (Sm) is introduced. This 
metric is calculated by averaging the square difference between two consecutive rows 
over the nodule candidate subimage. 

10 

J n m p 
Sm = S Z Z ( ^^^'J>^^ - ^ ^'^)>' 

nmp i J jfc 

The metric Sm is defined with respect to a global coordinate system as 
follows: 

15 the X axis runs between a patient's shoulders; 

the Y axis runs between a patient's chest and back; and 

the Z axis runs along the length of a patient's body, 
where i, j, and k are indices for X, Y and Z coordinates respectively and n,m and p are 
the dimensions of the subimage in X,Y and Z directions respectively. 

20 

A streaking artifact filter is applied for nodule candidate sub-images with a 
metric value greater than an empirically selected threshold. Tsar- The preferred range 
for Tsar is between 55,000 and 65,000 HU^. The preferred value of Tsar is about 
60,000 HU^, which corresponds to a pixel pair variation of about 245 HU. In the 
25 present invention, a median filter is implemented. Because of horizontal nature of the 
streaking artifact, vertical median filters are used. Preferably a vertical median filter 
of size 1x3 is selected. Figures 10 and 1 1 illustrate a nodule sub-image before and 
after streaking artifact removal. 

30 A multi-stage filtering procedure, illustrated in Figure 12, is used to reduce the 

number of false positive candidate regions. This approach is effective in minimizing 
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the computation time without sacrificing performance. The present invention 
preferably employs at least two filters. The first filter removes vessels and large 
vessel bifurcation points from the nodule candidate list. A large percentage of false 
positives are eliminated by the first filter. The nodule candidates list is then further 
refined by the second filter, which removes small bifurcation points. A detailed 
description of each filter is presented below. 

Nodule candidates generated in the h3^oAesis generation stage include 
nodules, blood vessels and bifurcation points which passed the AJR criterion. The 
objective of the first filtering stage is to remove blood vessels and large vessel 
bifurcation points firom the nodule candidate list. 

In this first filtering stage, the nodule candidate list is refined based on each 
nodule candidate's attachment surface area to other structures as shown in Figure 13. 
For this purpose, a fraction (F^) of the nodule surface that is attached to other solid 
structure is calculated by 

Fa = Sa/Sn (7) 

where Sn denotes the surface area of the nodule candidate and Sa denotes the 
attachment surface area to all other structures. 

Nodule candidates can be categorized into two groups based on the value 
of Ffl. The first group consists of vessels and large bifurcation points which have 
high values of Fa. The second group consists of nodules and very small 
bifurcation points. These nodule candidates have values of Fa smaller than the 
first group. A threshold can be used to eliminate the first group i.e vessels and 
large bifurcation points, based on the value of Fa. 

To estimate Fa, we need to find the outer surface of the nodule. This is 

A 

achieved by a region growing algorithm starting from R mi and incrementing by a 
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spherical layer until the anticipated R mc is reached. At this point. Fa is computed. 

A 

R MC is obtained from a precomputed table of values for a given Rmi- 

The current implementation of the procedure is described for a 2D region as follows. 

A 

5 1 . J? is obtained from a precomputed table of values for a given Ruff. 

2. For simplicity and efficiency, a rectangular rather than spherical region 
growing procedure is used. This is illustrated for a 2D region in Figure 14. 
Starting with a square region, the size is incremented by one pixel on all sides 
10 during every iteration. This iteration is repeated until the size of the square 

A 

matches R mo The pixels at the edges of the square region that are connected 
to the labeled nodule region are used to estimate 5^. The labeled nodule 
region is used to estimate Sn. The extension to three dimension is 
straightforward except for the anisotropic sampling space in low resolution CT 
1 5 scans. The 3D region growing algorithm increments at a different rate in the 

axial direction than in the in-plane direction. This growth rate is inversely 
proportional to the image resolution. For example, if the ratio of axial to in- 
plane resolution in an image is 3 : 1 , then three iterations of region growing is 
performed in the in-plane direction before growing once in the axial direction. 

A ^ 

20 3, Fa is estimated from iSfl/5„. 

A 

4. If Fa ^ Tay the nodule candidate is discarded. Ta is an empirically selected 
threshold representing the ratio of vessel attachment area to the total surface 
area of the nodule. Ta is preferably selected to be in a range between about .10 
25 and .45, and most preferably about .24. 

After the first filter, the majority of false positives are small (2mm size) vessel 

bifurcation points. Small vessel bifurcation points occur at junctions of small vessels. 

During global intensity thresholding procedure in the hypothesis generation stage, 

30 small vessels are often eliminated leaving the bifurcation point as a compact shape 
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region with very small attachment (as shown in Figure 15). This is due to the partial 
voxels effect. The nodule attachment analysis technique used in the previous stage 
generally does not identify these situations since the vessels were not present. 



5 The second filter eliminates bifurcation points from the nodule candidate list. 

A hollow cube with a certain wall thickness is used as shown in Figure 17. The 

A 

cube has an inner side length determined from nodule candidate radius R mi 
multiplied by a scale factor K and thickness 62 empirically selected. The preferred 

range of K is between 1 .5 and 2.5 and is preferably about 2. The thickness 82 
10 preferably has a minimum value of 1 and should be large enough to achieve noise 
insensitivity. The preferred range of 62 is between 2 and 6 voxels and is preferably 
4 voxels. The volume ratio, Ayr is then defined using the volume of the intersection 
Vni of the cube wall with the nodule candidate region and the nodule candidate 

A 

volume Vn as shown in equation 8. R mi is used to calculate F„. Bifurcation points 
1 5 have a higher Ayr value compared to small nodules. The threshold parameter used 
for this measure is TVv The preferred range for TVv is from about -774 to -674 HU, 
and is most preferably selected to be about -724 HU. 
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Ayr= VJVn (8) 
B. PULMONARY NODULE REGISTRATION 



Two high-resolution CT images are needed in order to measure the volume 
change of a nodule. The nodule is segmented from both CT images, and the percent 

25 volume change and the doubling time are calculated. With current methods, the 

segmentation of the nodule in the first image is independent from the segmentation of 
the nodule in the second image. The consistency of the nodule segmentation between 
the two images is improved by comparing the segmentation from the first image with 
the segmentation from the second image. 

30 A region-of-interest (ROI) is selected for the nodule in both of the CT images. 

The adaptive threshold is selected for each ROI using the histogram of the region, and 

both ROIs are then resampled to isotropic space. 
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To facilitate a meaningfiil comparison, the second ROI is registered to the first 
ROI. A rigid-body transformation model is assumed, meaning that in general the 
structures in the ROIs are confined to only translation and rotation. This is a valid 
assumption because the registration will be performed on a small focused area in the 
5 lungs. Thus, it is not expected that the region will change dramatically with the 
different levels of inspiration that may occur during the two CT scans. Registering 
the second nodule to the first nodule will change the location and orientation of the 
second nodule while preserving its volume and shape. 

10 After rigid-body registration, the nodule from the first image and the nodule 

from the second image will have the same orientation and position. The nodule is 
segmented firom first image and firom the registered second image using gray-level 
thresholding. The attached vessels are removed using iterative morphological 
filtering, and any pleural attachments axe removed using an iterative clipping plane. 

15 

The nodule segmentation of the two ROIs produces two binary images 
representing a nodule at two different times. The nodule firom the second image has 
been registered with the nodule in the first image. By looking at the corresponding 
pixels between the segmented nodules and the threshoddnd images, it is possible to 
20 label the pixels in the segmented nodules as nodule repeat pixel, nodule growth 

pixels, nodule atrophy pixels, or missegmented pixels. The nodule segmentations are 
adjusted by removing the missegmented pixels. This improves the consistency of the 
segmentations of the nodule in two different times, thus improving the accuracy of the 
volume measurements. 
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The algorithm for the registered segmentation procedure is as follows: 

1 . Select a Region-Of-Interest for the nodule in the two images. 

2. Adaptive thresholding selection. 

3. Isotropic resampling. 
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4. Register the second nodule to the first nodule using the rigid-body 
transformation. 

5. Perform segmentation on both of the nodules: 

(a) Gray-level thresholding using adaptive threshold. 

(b) Vessel Segmentation using iterative morphological filtering. 

(c) Pleural Wall Segmentation using a clipping plane. 

6. Rule-Based adjustment of the two nodule segmentations. 

The details of each step are explained below. 

A typical 3D high-resolution CT scan may contain anjwhere firom 5 to 30 
image slices. Each image slice has a resolution of 512 x 512 pixels. A small cubic 
region-of-interest with a size of about three times the diameter of the nodule is 
selected. Selecting a small ROI will reduce the amount of computation needed to 
register and segment the nodule. 

The segmentation of nodule tissue firom lung parenchyma can be achieved 
by either a fixed gray-level threshold or by an adaptive threshold. The threshold is 
selected by choosing a value that best separates the nodule (foreground) firom the 
lung parenchyma (background), which in the normal case is bimodal intensity 
histogram (see Figure 18). This separation may be achieved by selection of a 
threshold that falls at the midpoint between the peaks of each mode of the 
histogram. 

The fixed threshold is selected using the mean values of lung parenchyma 
and solid nodule tissue compiled over several cases. The mean value of lung 
parenchyma and solid nodule tissue was determined to be -880 HU and -27 HU, 
respectively. Thus, the fixed threshold was calculated to be -453 HU. 
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However, a fixed threshold does not take into account the change in lung 
parenchyma density due to inspiration of the lungs. A study of CT scans by the 
inventors determined that the density of the parenchyma around the nodule changed 
on average 9.7 HU ± stddev 7.0 HU (maximum of 21 HU) between repeat scans of 
the same patient. Furthermore, it was observed that the lung parenchyma density 
increases towards the back of the lungs when the patient is on his back because of 
the accumulation of blood due to gravity. In addition, a fixed threshold is not robust 
to changes in attenuation values due to differences in X-ray dosage and other CT 
scan parameters. 



Preferably an adaptive threshold is used to improve over the limitations of 
the fixed threshold. The present invention adaptively selects an optimal threshold 
for each nodule region-of-interest. Given a nodule region-of-interest, a histogram is 
calculated between -1024 HU and 476 HU using a bin size of 1 HU. The histogram 

15 may be noisy because of the small bin size. Accordingly, the histogram is 

smoothed by filtering with a Gaussian with a standard deviation of 25 HU. The 
peak values of parenchyma and solid nodule are expected to be relatively close to 
the mean values calculated over several cases. Thus, the peak parenchyma value 
calculated by searching over a range of 200 HU from the ideal value of -880 HU. 

20 Similarly, the peak solid nodule value is determined by searching over ±200 HU 
from the ideal value of -27 HU. Finally, the adaptive threshold is calculated as the 
midpoint between the two peak values. The adaptive threshold selection algorithm 
is as follows: 



25 1 . Given a nodule region-of-interest: 

2. Calculate the histogram, H{x), between -1024 HU and 476 HU with a 
bin size of 1 . 

3. Filter H{x) with a Gaussian with stddev of 25 HU. 

4. Peak Parenchyma: Vp = max-\024<x<'6^Q H{x) 
30 5 . Peak Nodule : Vn = max-227<;x<i73 H{x) 

6. Adaptive Threshold: thay = (v^ + v„)/2 
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Figure 18 shows an example histogram of a 12mm nodule. The parenchyma 
peak and solid tissue peak were determined to be -848 HU and -16 HU respectively. 
The adaptive threshold was then calculated as -432 HU, as opposed to -453 HU for 
the fixed threshold. The histogram around the solid-tissue peak is seen in Figure 21. 
5 Figures 20 and 21 shows the histogram of an ideal 12mm nodule. The histogram of 
the real nodule has more voxels in the transition region between parenchyma and 
solid tissue because it contains partial voxels from small vessels and other objects. 

For small nodules less than 5mm in size, the solid-tissue distribution in the 
histogram is usually obscured by partial voxels from small vessels and noise voxels. 

10 In this case, the solid-tissue peak value cannot be found through the histogram. The 
adaptive thresholding algorithm is modified such that the mean solid-tissue value is 
fixed at the model value of -27 HU. The adaptive threshold is then calculated as the 
midpoint between the peak parenchyma value and -27 HU. A histogram of a 5mm 
nodule appears in Figure 22. The mean solid-tissue value cannot be found because 

1 5 there is no prominent peak on the right side of the histogram. 

After the adaptive threshold is determined, the ROI is then resampled into 
0.25mm isoptropic space using trilinear interpolation. 

20 The second ROI is registered to the first ROI using a rigid-body 

transformation. The rigid-body transfor-mation is a function of six parameters: (tx, 
tyy tz, rx, ry, r^, or translation in x, translation in y, translation in z, rotation about the 
X-axis, rotation about the y-axis, and rotation about he z-axis, respectively. The 
rigid-body transformation is a mapping of a point v in 3-d space to a point v' in 

25 transformed space defined by the following equation: 



(1) 
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where Rx, Ry, and Rz are the rotation matrices defined as: 
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R. = 



R, = 
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0 cos(r,) -sin(r^) 
0 sin(r,) cos(r,) 
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(2) 



(3) 



(4) 
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The second ROI is registered to the first ROI by conducting a search over the 6 
transformation parameters: The search is designed to minimize a similarity metric 
between the transformed second ROI and the first ROI, thus aligning the two ROI's 
using the transformation parameters. The mean-squared-difference (MSD) is used 
as the similarity metric and is defined as 



^■5^=T7 ZZZ iROI^ij,k)-ROhiijM 



(5) 



i.J,k 
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where N is the number of pixels in ROI\. With the MSD metric, perfectly registered 
images will produce a metric of zero. The goal of registration is to align the nodule 
in the firet ROI with the nodule in the second ROI. The MSD metric can be 
modified by weighting it with an appropriately sized gaussian. The gaussian is 
positioned over the center of the nodule and has a standard deviation equal to the 
radius of the nodule. The gaussian-weighted MSD is defined as: 



exp 



iROUijJC)-ROHijJk)f (6) 
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where (gxigy^gz) is the center of the gaussian and gr is the standard deviation of the 
gaussian. Using the gaussian-weighted MSD will force the registration algorithm to 
focus more on aligning the nodule than any periphery structures (e.g. vessels or 
pleural wall). However, some of the periphery structures are needed for a good 
5 registration, thus the size of the gaussian is selected to include all of the nodule and 
some of its periphery. 

The similarity metric is minimized using Powell's method, a multi- 
dimensional direction-set search algorithm. W. Press, editor. Numerical Recipes in 

10 C. Cambridge University Press, second edition, 1992. The initial translation 

parameter is set to the offset of the center of the first nodule to the center of the 
second nodule. The initial rotation parameters are set to (0,0,0). The search is 
terminated when either all the transformation parameters change within some 
epsilon, or the similarity metric does not change by more than a tolerance value. 

15 After the registration is complete, the second ROI will be the same size as the first 
ROI. 

An example of rigid-body registration on two nodules is shown in Figure 26. 
The first ROI and the second ROI are shown in Figures 24 and 25, respectively. 

20 The registration algorithm aligned the second ROI with the first ROI using 

translation parameters of (-13.40,16.61, -0.88) pixels and rotation parameters of 
(0.79, -3.52,-7.20) degrees. The resulting registered image is seen in Figure 26. 
Figure 27 shows the difference image between the first ROI and the second ROI. 
The difference image shows that the vessels at the top of the image and on the right 

25 of the image are not aligned properly. The circular white ring is the growth of the 
nodule between the two scan times. Figure 28 shows the difference image between 
the first ROI and the registered second ROI. The vessels are less visible in the 
second difference image than in the first difference image, meaning that with the 
registered image the vessels and nodule are aligned better than without registration. 

30 

A binary image of the solid nodule is created by thresholding the image with 

the adaptive threshold determined previously. Any attached vessels are removed by 
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using iterative morphological filtering. A morphological opening, followed by 
masking with the original binary image, is performed using a decreasing sphere size. 
The initial diameter of the sphere is 3/4 of the nodule radius, resulting in the 
removal of all vessels that are smaller than 3/4 of the radius. 

If the nodule is juxtapleural, then the pleural surface is segmented from the 
nodule using a clipping plane approach. Starting at a point inside the nodule, a 
clipping plane is moved towards the pleural wall. When the relative change in the 
size of the nodule is over a threshold, the plane is reorientated to minimize the size 
of the nodule. Finally, the algorithm stops when the reorientation procedure does 
not produce any significant change in the size of the nodule. The juxtapleural 
nodule is then segmented by using the clipping plane from the previous iteration. 

Figure 29 shows an example of the solid nodule segmentation. A montage of 
the original gray-scale nodule is shown on in Figure 29, while the segmented nodule 
is shown in Figure 30. The segmented nodule is shaded in red. The pleural surface 
and vessels are shaded in green, while the background is shaded in gray. The solid 
nodule segmentation has detached the vessel and pleural wall from the nodule. 

Given the two segmented nodules, it is possible to adjust the nodule 
segmentations by comparing the segmented nodules and the thresholded regions. 
Let iSi be the segmented nodule from the first image and T\ be the thresholded first 
image before vessel or pleural surface removal. Likewise, let iS'2 be the segmented 
nodule from the second image and T2 be the thresholded second image. 

A rule-based system is used to mark active pixels in the segmented nodule S\ 
as repeat nodule, nodule atrophy, or nodule missegmentation. If an active pixel in 
the first segmented nodule corresponds to an active pixel in the second segmented 
nodule, then that pixel is a repeat nodule pixel because it is present in both 
segmented nodules. If an active pixel in the first segmented nodule corresponds to 
an inactive pixel in the second segmented nodule and an inactive pixel in the second 
thresholded image, then that pixel is a nodule atrophy pixel because it is present in 
the first nodule, but not present in the second image. Finally, if an active pixel in 

37 



the first segmented nodule corresponds to an inactive pixel in the second segmented 
nodule and an active pixel in the second thresholded region, then that pixel is a 
missegmented nodule pixel because it is a non-nodule object (vessel or pleural wall) 
in the second region. 

An active pixel in the segmented nodule 52 can be marked as repeat nodule 
pixel, nodule growth pixel, or nodule missegmented pixel by using a similar set of 
rules. A summary of the rules for marking a region in 5i or 52 are given below. 

For a foreground pixel in the first segmented nodule, S\: 

• (A) If the corresponding pixel in 52 and the corresponding pixel in T2 are both 
foreground, then the pixel in 5i is a repeated nodule pixel from the first 
region-of-interest (ROIi) to the transformed second region-of-interest (ROl2t). 

• (B) If the corresponding pixel in 52 is background and the corresponding pixel 
in T2 is background, then the pixel in S\ is nodule atrophy pixel. 

• (C) If the corresponding pixel in 52 is background and the corresponding pixel 
in T2 is foreground, then the pixel in S\ is a missegmented pixel in the first 
region-of-interest (ROIi). 

For a foreground pixel in the second segmented nodule, 52: 

• (D) If the corresponding pixel in S\ and the corresponding pixel in T\ are both 
foreground, then the pixel in 52 is a repeated nodule pixel from the first 
region-of-interest (ROIi) to the transformed second region-of-interest (ROtt). 

• (E) If the corresponding pixel in S\ is background and the corresponding pixel 
in Ti is background,, then the pixel in 52 is nodule growth pixel. 

• (F) If the corresponding pixel in S\ is background and the corresponding pixel 
in Ti is foreground, then the pixel in 52 is a missegmented pixel in the 
transformed second region-of-interest (ROht). 

Figures 31 through 38 illustrate an example of the segmentation adjustment on 
a registered vascularized nodule. The thresholded nodule at two different times is 
shown in Figures 31 and 32, and the segmented nodules are shown in Figures 33 and 
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34. In the first segmentation part of the top attached vessel has been missegmented as 
part of the nodule, and in the second segmentation, part of the bottom attached vessel 
has been missegmented as part of the second nodule. By using the rules described 
above, the regions of each segmentation are marked in Figures 35 and 36 as (A) 
5 nodule in time 1 , (B) nodule atrophy, (C) nodule missegmentation in time 1, (D) 
nodule in time 2, (E) nodule growth, or (F) nodule missegmentation in time 2. 
Finally, the segmentations of both nodules are adjusted in Figures 37 and 38 by 
removing the missegmented regions (C) and (F). The volumes of the nodules in the 
adjusted segmentations are more accurate than the volumes in the original 
10 segmentations because the vessel attachments have been removed. This leads to a 
more accurate determination of doubling time or percent volume change. 
Furthermore, regions of growth and atrophy of the nodule can be examined over time. 

Thus, while there have been described what are presently believed to be the 
1 5 preferred embodiments of the invention, those skilled in the art will realize that 

changes and modifications may be made thereto without departing firom ttie spirit of 
the invention, and is intended to claim all such changes and modifications as fall 
within the true scope of the invention. 
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